arXiv: 1503.03705v4 [q-fin.CP] 27 Jan 2016 


A hybrid tree/finite-difference approach for 
Heston-Hull-White type models 

Maya Brian:* 

Lucia Caramellino^ 

Antonino Zanette^ 

January 29, 2016 


Abstract 

We study a hybrid tree/finite-difference method which permits to obtain efficient and accurate 
European and American option prices in the Heston Hull-White and Heston Hull-White2d mod¬ 
els. Moreover, as a by-product, we provide a new simulation scheme to be used for Monte Carlo 
evaluations. Numerical results show the reliability and the efficiency of the proposed methods. 
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1 Introduction 

In this paper we consider the Heston-Hull-White model, which is a joint evolution for the equity value 
with a Heston-like stochastic volatility and a generalized Hull-White stochastic interest rate model 
which is consistent with the term structure of the interest rates. We consider a further situation where 
the dividend rate is stochastic, a case which is called here the “Heston Hull-White2d model” and can 
be of interest in the multi-currency (the dividend rate being interpreted as a further interest rate). We 
concern the problem of numerically pricing European and American options written on these models. 

At the present time, the literature on this subject is quite poor and includes Eourier-Cosine meth¬ 
ods, semi-closed approximations and finite-difference methods to price vanilla options. In [15], Grzelak 
and Oosterlee introduce two approximations of the non-affine models. The Eourier-Cosine method is 
then used on this approximating affine model. The authors remark that for accurate modeling of 
hybrid derivatives it is necessary to be able to describe a non-zero correlation between the processes 
driving the equity and the interest rate. This is possible in the approximations presented in their 
paper but only using approximated affine models. Haentjens and in’t Hout propose in [13] a finite- 
difference Alternating Direction Implicit (ADI) scheme for pricing European options solving the orig¬ 
inal three-dimensional Heston-Hull-White Partial Differential Equation (hereafter PDE). The Heston 
Hull-White2d model is treated using semi-closed approximations in the Foreign Exchange model HE]. 

In this paper, we generalize the hybrid tree/finite-difference approach that has been introduced 
for the Heston model in the paper [6]. In practice, this means to write down an algorithm to price 
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European and American options by means of a backward induction that works following a finite- 
difference PDE method in the direction of the share process and following a recombining binomial 
tree method in the direction of the other random sources (volatility, interest rate and possibly dividend 
rate). 

It is well known that there is a link between tree methods and finite-difference methods. The 
most remarkable benefits in using recombining binomial trees (let us stress the terms “recombining” 
and “binomial”: just two possible recombining jumps at each time-step for each component) are the 
simplicity of the implementation, the low computational costs and the efficiency of the output numer¬ 
ical results. In dimension 1, one can always build a recombining binomial tree (see e.g. Nelson and 
Ramaswamy |20]) but this is not the case in multidimensional problems. For example, in the standard 
(dimension 2) Heston model it is not possible to write down a recombining binomial approximating 
tree - roughly speaking, this follows from the fact that it is not possible to produce a function of the 
Heston components such that the associated Stochastic Differential Equation (SDE) has a diagonal 
diffusion coefficient. A binomial tree approximation for the standard Heston model has been proposed 
by Vellekoop and Nieuwenhuis in |22j but it is far from being recombining and, as shown in [6], it is 
problematic when the Feller condition is not satisfied. Finally, an approximation of the price coming 
from a numerical treatment of the (multidimensional) PDE can be very expensive, mainly to handle 
the 4-dimensional Heston-Hull-White2d model. 

So, the idea underlying the approach developed in this paper is in some sense very simple: we apply 
the most efficient and easy to implement method whenever we can do it. In fact, wherever an efficient 
recombining binomial tree scheme can be settled (volatility, interest rate and possibly dividend rate), 
we use it. And where it cannot (share process), we use a standard (and efficient, being in dimension 1) 
numerical PDE approach. Hence we avoid to work with expensive (because non recombining and/or 
binomial) trees or with PDFs in high dimension. Moreover, for the Cox-Ingersoll-Ross (hereafter CIR) 
volatility component, we apply the recombining binomial tree method firstly introduced in [1], which 
theoretically converges and efficiently works in practice also when the Feller condition fails. 

The description of the approximating processes coming from our hybrid tree/hnite-difference ap¬ 
proach, suggests a simple way to simulate paths from the Heston-Hull-White models. Therefore, we 
propose here also a new Monte Carlo algorithm for pricing options which seems to be a real alternative 
to the Monte Carlo method that makes use of the efficient simulations provided by Alfonsi [T]. 

Our approaches allow one to price options in the original Heston-Hull-White processes with non¬ 
zero correlations. Here, we consider the case of a non null correlation between the equity and the 
interest rate process, as well as between the equity and the stochastic volatility. Moreover, in the 
Heston-Hull-White2d model, we allow the dividend rate to be stochastic and correlated to the equity 
process. But it is worth noticing that other sets of correlations can surely be selected. 

The paper is organized as follows. In Section we introduce the Heston-Hull-White model. Then 
in Section we construct a recombining binomial tree approximation for the pair given by the volatility 
and the interest rate process. Section refers to the approximation of functions of the underlying 
asset price process by means of PDE arguments. In Section we describe the hybrid tree/hnite- 
difference scheme for the computation of American options. In Section we see how to generalize the 
previous procedure in order to handle the Heston-Hull-White2d process. In Sectionwe show that our 
arguments can be used also to set-up simulations, to be applied to construct Monte Carlo algorithms. 
Finally, numerical results and comparisons with other existing methods are given in Section showing 
the efficiency of the proposed methods in terms of the results and of the computational time costs. 
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2 The Heston-Hull-White model 


The Heston Hull-White model concerns with cases where the volatility V and the interest rate r are 
assumed to be stochastic. The dynamics under the risk neutral measure of the share price S and the 
volatility process V are governed by the stochastic differential equation system 

^ = {rt- r])dt + \/Vt dZt, 

Jt 

dVt = KviOv - Vt)dt + avVVtdWl, 
drt = Kr{9r{t) — rt)dt -|- ardW^, 

with initial data 5o > 0, Vq > 0 and ro > 0, where Z, and are suitable and possibly correlated 
Brownian motions. Recall that 14 is a CIR process whereas rt is a generalized Ornstein-Uhlenbeck 
(hereafter OU) process: here Or is not constant but it is a deterministic function which is completely 
determined by the market values of the zero-coupon bonds (see [ZD- 

Let us fix the correlations among the Brownian motions. As observed in [15], the important 
correlations are between the pairs (S, V) and {S,r). So, we assume that W = (W^, W’^) is a standard 
Brownian motion in and Z is a Brownian motion in M which is correlated both with and W'^: 

d{Z, Wi)t = Pi dt and d{Z, W 2 )t = P 2 dt. 

By passing to the logarithm Y = In S' in the first component and taking into account the above 
mentioned correlations, we reduce to the dynamics 

dYt = {rt-v- ^Vt)dt + VVt {pidWl + p 2 dW^ + psdWj ^), Tq = In Sq G M, 

dVt = Kv{Ov -Vt)dt + avVVt dW/, Rq > 0, 

drt = Kr{9r{t) — rt)dt -|- ardWt, vq > 0, 

where W = iW^, IT^, W^) is a standard Brownian motion in and the correlation parameter p^ is 
given by 

P3 = pI- pI with p\ + pI< 1. 

As already done in m, the process r can be written in the following way: 


rt = OrXt -I- pt 


where 


Xt — Kr 


rt 

Xgds+Wf and pt = + Kr / 9r{s)e~^"'d-^^ds. 

Jo 


So, we can consider the triple (Y, V,X), whose dynamics is given by 


dYt = pY{Vt, Xt, t)dt + VR {pidWl + p2dWi + p^dWf ), To = In So G M, 
dVt = pv{Vt)dt + avVVt dWl, Ro > 0, 
dXt = px{Xt)dt + dWf, Xo = 0, 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 
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where 


firiv, X, t) = arx + ipt-ri-]^v. 

(2.4) 

fiv{v) = Kv{0v - v). 

(2.5) 

flx{x) = —KrX. 

(2.6) 


The purpose of this paper is to efficiently approximate the process {Y, V,X) in order to numerically 
compute the price of options written on the share process S. 


3 The recombining binomial tree for the pair X and V 

First of all, we consider an approximation for the pair {V,X) on the time-interval [0,r] by means of 
a 2-dimensional computationally simple tree, that is by means of a Markov chain that runs over a 2- 
dimensional recombining bivariate lattice (recombining binomial tree). In the usual case, as in the Cox- 
Ross-Rubinstein tree m, at each time step the process can jump either on the nearest up-node or on 
the nearest down-node. Here, we consider the possibility of “multiple jumps” as introduced in Nelson 
and Ramaswamy [20]. Roughly speaking, the process can again jump upward or downward but the 
up/down jump nodes might not be the nearest ones: they are dehned as the up/down positions at the 
next time-step whose associated transition probabilities better interpolate the theoretical expectation 
of the transition. As discussed in Nelson and Ramaswamy [20) . this is the best way to construct 
an efficient tree for the approximation of one-dimensional diffusion processes, especially when the 
diffusion coefficient is not constant. Figure shows an example of possible “multiple jumps” for the 
trees that approximate our processes X and V, that we are going to describe. 

In this section, we consider a discretization of the time-interval [0, T] in N subintervals [nh, (n+l)h], 
n = 0,1,..., N, with h = T/N. 


3.1 The tree for X 

The construction of the recombining binomial tree for the process X is quite standard, because here 
the diffusion coefficient is constant. For n = 0,1,..., A^, consider the lattice for the process X 


xll = {xn,j}j=o,i,...,n with Xn,j = (2j - n)Vh 


(3.1) 


(notice that xo,o = 0 = Xq). For each fixed Xnj G we define the “up” and “down” jump by means 
of Ju(n,j) and ja{n,j) defined by 

ju{n,j) = min{j* : j + l< 3* <n+l and Xn,j + tixixn,j)h < Xn+ij*}, (3.2) 

JdinJ) = max{j* : 0 < j* < j and Xn,j + fix{xn,j)h > Xn+i,j*}, (3.3) 


fix being the drift of the process X, see (2.6). As usual, one sets ju{n,j) = n-|- 1 if {j* : J + 1 < j* < 
n + 1 and Xn,j + fixixn,j)h < = 0 and j^{n,j) = 0 if {j* : 0 < j* < j and Xn,j + fix{xn,j)h > 

Xn+i,j*} = 0- Note that the up/down jumps in ( |3.2| )-( [3^ might not be the nearest up/down positions 
in the lattice at time n -|- 1. An example is given in Figure [^left, where the lattice X^ is drawn and 
some possible instances of Xnj, x^_^_l and shown to exhibit as the tree can be 

visited. 
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The transition probabilities are defined in order to better interpolate the expected local transition: 
starting from the node (n,j), the probability that the process jumps to j^(n, j) and j^{n,j) at time- 
step n + 1 are set as 


Pu' {n,j) = Q\/ --- ' ^ ^ A 1 and pj’ (n,j) = 1 


(3.4) 


respectively. This gives rise to a Markov chain (^^)n=o,...,Af that weakly converges, as /i —>■ 0, to the 
diffusion process (-^t)t6[o,T] turns out to be a robust tree approximation for the OU process X. 


3.2 The tree for V 


For the CIR volatility process V, we consider a recombining binomial tree procedure that again 
follows the “multiple jumps” approach. In this case the recombining lattice is built by means of the 
transformation 

m) = — Vvt. 

cry 

This transformation is particularly important because /(Vj) turns out to be a diffusion process with 
unit diffusion coefficient, and this fact is useful in order to construct a recombining lattice. Many 
authors (see e.g. m or [23]) propose tree algorithms for Vj by working on the transformed process 
f{Vt). The unpleasant fact is that now the drift of /(Vj) is very bad and is such that the approximating 
process converges only when the Feller condition holds: 2kv0v ^ cry- In order to overcome this fact, 
we use the approach in [Ij, that, roughly speaking, works as follows: the tree structure is built by using 
again / (see next ( |3.5[ )) but the possible jumps and the transition probabilities are set on the dynamics 
of the original (and not transformed) CIR process V) (see next (3.6)-(3.7) and (3.8)). The main fact 
is that now the weak convergence on the path space is achieved for every values of ny^Oy^ay > 0, 
so the Feller condition is not required. Details and comparisons with other tree existing methods to 
approximate the CIR process are given in |3]. 

For n = 0,1,..., A^, consider the lattice 


k=0,l,...,n with 


Vo + ^{2k-n 




\/Vo+ % {2k—n)\/Ji>0 


(3.5) 


(notice that uo,o = Cq). For each fixed Vn,k £ V^, we define the “up” and “down” jump by means of 

k'^{n, k) = : k + 1 < k* < n + 1 and Vn,k + Pv{vn,k)h < Vn+i,k*}, (3.6) 

k^{n, k) = max{/c* : 0 < k* < k and Vn,k + Pvivn,k)h > v^+i^k*} (3-7) 


where the drift fty of V is dehned in (2.6) and with the understanding k^{n,k) = re -|- 1 if {k* : 
k + l<k*<n + l and Vn,k + py{vn,k)h < Vn+i,k*} = 0 and k^{n,k) = 0 if {k* : 0 < k* < 
k and ^ -|- ^y{vn,k)h > Vn+i,k*} = 0- By construction, the up/down jumps in (|3.6|)-(3.7) might not 


be the nearest up/down positions in the lattice at time re + 1. In Figure [fright we show an example 
of the lattice together with some possible instances of the triple k^(n j)^ '^n+i,k>^(n,j))■ 

The transition probabilities are defined in order to better interpolate the expected local transition: 
starting from the node (re, k) the probability that the process jumps to k^{n, k) and k^{n, k) at time- 
step re -b 1 are set as 


Vh/ pv{Vn,k)h + Vn,k ^n+l,fcRn,fc) . , yh, ^ yh/ 

p)(’"(re, fe) = 0 V- - - -/\1 and p^’^{n,k) = I — p^’^{n,k) 

'^n+l,k!i^{n,k) ~ '^n+l,k^{n,k) 


(3.8) 
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respectively. This gives rise to a Markov chain that weakly converges, as /i —)• 0, to the 

diffusion process (Vt)tG[o,T] and turns out to be a robust tree approximation for the CIR process V - 
details are given in [3]. 
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Figure 1: Example of a tree for X on the left and of a tree for V on the right. 


3.3 The tree for the pair (V,X) 

The tree procedure for the pair {V,X) is set by joining the trees built for V and for X. Namely, for 
n = 0,1,..., N, consider the lattice 




(3.9) 


Starting from the node {n,k,j), which corresponds to the position {vn,k,Xn,j) G x X^, we define 
the four possible jump by setting the four nodes at time n + 1 following the definitions (3.2)-(3.3) and 


(3.6)-(3.7): 


{n + l,k^{n,k),jii{n,j)) with probability Puuin,k,j) = pu'^{n,k)pi’'"{n,j), 

(n +l,/c()(n,/c),j^(n,j)) with probability p'ld{n,k,j) = Pu’^{n,k)p^'^{n,j), 

{n + l,k^{n,k),j^{n,j)) with probability p^^{n,k,j) = p'^’^{n,k)pu'^{n,j), 

(n +l,A:^(n,A;),j^(n,j)) with probability p^^{n,k,j) = p^'^{n, k)p^’^{n, j). 


x,h, 


(3.10) 


where the above probabilities pu'^{n,k), p^'^{n,k), Pu'^{n,j) and p^'^{n,j) are defined in (3.8) and 
(3.4) respectively. The above factorization is due to the orthogonality of the noises driving the two 
processes. As a quite immediate consequence of standard results (see e.g. the techniques in m) , one 
gets the following; the associated bivariate Markov chain (tv, ^^)n=o,...,Af weakly converges to the 
diffusion pair {Vt, solution to 


dVt = fiv{Vt)dt + av^/VtdWl, Eo > 0, 
dXt = —KrXt dt + CJr dW^, Xq = 0. 


Remark 3.1 In the case one is interested in introducing a correlation between the noises and W'^ 
driving the process V and X respectively, the joint tree can he constructed on the same lattice but the 
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jump probabilities are no more of a product-type: the transition probabilities 

p^^{n,k, j) can be computed by matching (at the first order in h) the conditional mean 
and the conditional covariance between the continuous and the discrete processes of V and X. More 
precisely, for both components the conditional mean is matched by construction (this is actually the 
main consequence of the definition of the multiple jumps). As for the conditional covariance, assuming 
that d{W^,W‘^)t = adt, with |a| < 1, then one has d{V,X)t = aavV^dt. Therefore, the matching 
conditions lead to solving the following system: 

' Puu{n,k,j)-^p^fin,k,j) =pu’^{n,k) 

Puuin,k,j) +p’j^in,k,j) =p^'^{n,j) 

' pLin,k,j)+p’(fin,k,j) + Pduin,k,j) + p'j^{n,k,j) = 1 
muuin, k,j)p^^{n, k,j) + k, j)p^a{n, k, j)+ 

+m'ju{n,k,j)p'ju{n,k,j)+m’j^{n,k,j)p’jfin,k,j) = aav^JUnfkh 

where 

i'^n+l,k!^{n,k) (^n+1 j5(n,j) 

~ ~ — Xn,j), 

~ ('^n+l,A:J(n,fc) “ (®n+lji(nj) “ Xn,j), 

j) ~ (^n+l,fcJ(n,A:) “ jJ(nj) “ ®nj)- 

This is done in in a different context but the proof of the weak convergence on the path space is 
analogous - this can be done by standard arguments, as in m or Uf. 


4 Approximating the y-component: the finite-difference approach 

We go now back to ( |2.3[ ), that is 

dYt = fiYiVt, Xt, t)dt + ^/y {pidWl + P 2 dwi + p^dWf ), Yo = In 5o, 
dVt = pv{Vt)dt + avVVt dWl, Fo > 0, 
dXt = px {Xt)dt + dttf, Xo = 0, 


where py-, Pv and px are given in (2.4), (2.5) and (2.6) respectively. By isolating y/VtdWf in the 
second line and dWf in the third one, we obtain 


dYt = ^dVt + p2VVtdXt + p{Vt, Xt, t)dt + P3^/VtdW^ 

(XV 


(4.1) 


with 

p{v,X,t) = pY{v,X,t) - -^pv{v) - P 2 \/v px{x) 

= (TrX + Tt - P - \ v - ^ Kv{0v - u) + P 2 nrXy/v. 

What we are going to do is mainly based on the fact that the noise is independent of the processes 
V and X. 
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4.1 The approximating scheme for the triple (Y,V,X) 

We consider an approximating process for Y turning out by freezing the coefficients in (4.1): we 
define Yq = Yq and for t G [n/i, (n + l)/i] with n = 0,1,..., iV — 1 we set 

yl" =Yl + —{Vt- Vnh) + P 2 VKh{Xt - X^h) + p{Vnh, Xnh, nh){t - uh) + 
av 

We consider now the approximating tree (IV) ^n)ne{o,...,Af} we call (1^^,the associ¬ 
ated time-continuous approximating process for the pair (F, X), that is 

We then assume that the noise driving the pair ^/^)te[o,T] is independent of the Brownian motion 
and we insert this discretization for (V, X) in the discretization scheme for Y. So, we obtain our 
final approximating process Y^^ by setting Y^ = Yq and for t G [nh, (n -|- l)h] with n = 0,1,..., N — 1 
then 


yt=yl+—{Vt-vlj 

av 


^hi^t-XL)+PiXL,Vl,nh){t-nh)+p3^/v^f^{W^ (4.3) 


Notice that if we set 


h_^h Pinrh -/vMx!^-Xnh), te[nh,{n + l)h] 


Zn ^ V 

av 


then we have 


dZ^ = nh)dt + dWl t G {nh, (n + l)h], 


yh _ -vh 

^nh ^nh 


(4.4) 


(4.5) 


that is solves a SDE with constant coefficients and at time nh it starts from 1^^. Take now a 
function /: we are interested in approximating 

®'(/(h^(n+l)/i) I Ynh ViYnh '^jX^h ^)- 


By using our scheme and the process Z^ in (4.4), we approximate it with the expectation done on the 
approximating processes, that is 

E(/(T(t+l)h) I Ynh = y. Vl = V, = x) 


= Hf{Z^n+l)h + ^(^(n+l)h - Vnh) + P2^JVj^hiXfn+l)h " I ^nh = V,X^h = ^) • 

Since {V^,X^) is independent of the Brownian noise driving Z^ in (4.4), we can write 


nf{ytn+l)h) I Ynh = y, VHt, = V, Xii^ = X) 


= - ^) + P2\Al(^(n+l)h -x)]y,v,-. 


vi = v.x':^^ = x]. 


in which 


4/p(e; y, V, x) = nf{Z\n+i)h + 0 I Z^h = y. K'). = X^h = x). 

Now, in order to compute the above quantity 'hp(^), consider a generic function g and set 

u{s,z;v,x) = E{g{Z^^^^^^) \ Z^ = z,V^ = v,X^ = x), s G [nh, {n + l)h]. 


(4.6) 

(4.7) 
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By (4.5) and the Feynman-Kac representation formula we can state that, for every fixed x G M and 


X > 0 , the function (s, z) i—)• ri(s, z] v, x) is the solution to 


dgU + /i(u, X, s)dzU + = 0, s G [nh, (n + l)h), z G M, 

u{{n + l)h,z;v,x) = g{z), 


(4.8) 


p being given in (4.2). In order to solve the PDE problem (4.8), we use a finite-difference approach 


4.2 Finite-differences 


At each time step n we numerically solve (4.8) at time s = nh by applying finite-difference techniques. 


We fix a grid on the y-axis = {Ui = Yq + with Jm = {-M, ..., M} and Ay = 


yi — yi-i. For fixed n, v >0 and x G M, we set = u{nh,yi;v,x) the discrete solution of (4.8) at time 
nh on the point yi of the grid Ym - for simplicity of notations, we do not stress in the dependence 
on V and x (from the coefficients of the PDE). 

The finite difference method we are going to set is inspired from the one developed in [H]. But a 
main difference arises: here, we do not distinguish anymore between the diffusion dominant or reaction 
dominant case and we propose to apply a full implicit finite-difference approximation in time. In fact, 
the discrete solution to problem (4.8) at time nh is computed in terms of the solution at time 


(n -|- l)h by using the following finite-difference scheme: 


- u. 


h 


— + p{v, X, nh) 


K+i - 


1 


i-l , ^ 2 




Ay2 


(4.9) 


Of course, (4.9) has to be coupled with suitable numerical boundary relations. We assume that the 


boundary values are defined by the following Neumann-type conditions: 


^-M-l — ^-M-l-1) '^M+1 — ^M-1- 


(4.10) 


Then, by applying the implicit finite-difference (4.9) coupled with the boundary conditions (4.10), 


we get the solution ..., by solving the following linear system 

Au^ = u^+\ 

where A = A{v,x) is the (2M -|- 1) x (2M -|- 1) tridiagonal real matrix given by 

\ 


(4.11) 


/ 1 + 2/3 -215 

—j3 + a 1 + 2,j3 —f5 — a 


A = 


(4.12) 


with 


a = 


—/3 + a 1 + 2j5 —f5 — a 
-2/3 1 + 2/3 y 


p{v,x,nh) and /3 = 777 —^ 


2Ay 


2Ay2 


(4.13) 


p being defined in (4.2). We stress that at each time step n, the quantities v and x are constant 


and known values (defined by the tree procedure for the pair {V,X)) and then a and /3 are constant 
parameters too. 
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One can easily see that the implicit scheme (4.9) is unconditionally stable. Moreover, by applying 


standard results (Theorem 2.1 in [ 8 ] e.g.), the matrix A is invertible for /3 7 ^ |a|. Therefore, setting 

n(u, x) = A~^(v, x), (4.14) 


the numerical solution to (4.8) on the grid Tm through the above discretization procedure is given by 

u(nh,yi;v,x) uf = ni£(v,x)g(ze), ieJM- (4.15) 

Remark 4.1 Other numerical boundary conditions can surely be selected, for example the two bound¬ 
ary values and may be a priori fixed by a known constant (this procedure typically appears 
in financial problems). 

4.3 The scheme on the T-component 

We can now come back to our original problem, that is the computation of the function ^f{f', y, v, x) 
in (4.7) allowing one to numerically compute the expectation in (4.6). 

We consider the approximating process {Y^, X^) as described in Section 4.1, This means that 

the pair (u, x) at time-step n is located on the lattice x X))-. v = Vn^k and x = Xn,j, for 0 < k, j < n. 
Then ( |4.15| ) gives the following approximation: for each yi G Ym, 

f yiXn,ki — j)f{y£YC), 

£&Jm 


Therefore, the expectation in (4.6) is computed on the approximating tree for {V,X) by means of the 
above approximation: 

®'(/(^(n-|-l)h) I ^nh ~ yii^nh ~ ''^n,ki^nh ~ ^n,j) — ^il{'£'n,ki Xn,j)Y n,k,j f {(-i a,b)p^^,{n,k, j) 

a,bG{d,u} I^Jm 

(4.16) 

where 

Yn,k,jfi(^,a,b) = f(^yi + ^{Vn+l,ka{n,k) - v) + P2y/viXn+l,j,,{n,j) “ 3 ^)) 


and the jump probabilities p^i,[n, k,j) are given in (3.10) (or in Remark 
between the noises driving V and X). 


3.1 


if a correlation is assumed 


Similar arguments can be used in order to compute the conditional expectation in the left hand 


side of (4.16) when the function / depends on the variables v and x also. Then one gets 


K+l)h.Y+l)k’ ) I Yk = », Vih = = Xnj) 

- E {vn,k,Xn,j)Tn,k,jf{^, a, h)p^^{n, k, j) 


(4.17) 


a,bg{d,n} I&Jm 


where 


T^n,k,jf{^,a,b) = 

f (fJi T ~ (.Xfi+l, ka{n,k) W,fc) T P2-\JVn,k jj ) > W-|-l,A:a(n,A:)) 


(4.18) 
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5 The algorithm for the pricing of American options 


The natural application of the hybrid tree/finite-difference approach arises in the pricing of American 
options. Consider an American option with maturity T and payoff function First of all, 

we consider the log-price process, so the obstacle will be given by 

= tG[0,T]. 

The price P{t, y, v, x) of such an American option is given by (recall the relation between the interest 
rate r and the process X: rt = OrXt + (ft, see (2.1)) 

P(t, 2 /,u,x)= sup 

rSTt.T ^ ^ 

where 7t,T denotes the set of all stopping times taking values on [t,T]. Hereafter, , X^’^) 

denotes the solution of the SDE (2.3) starting at {y,v,x) at time t. 

The price at time 0 of such an option is then approximated by a backward dynamic programming 
algorithm. Consider a discretization of the time interval [0, T] into N subintervals of length h = T /N: 
[0, T] = [J^~Q[nh, {n + l)h]. Then P(0, Tq) h); Aio) is numerically approximated through the quantity 
Ph{0, Yq, Vo,Xq) which is iteratively defined as follows: for (y, u, x) G M x M+ x M, 

J Ph{T,y,v,x) = T(y) and as n = iV - 1,..., 0 

\Ph{nK y, V, x) = max {T(y), + l)h, ^(t+i) j) }• 

From the financial point of view, this means to allow the exercise at the fixed dates nh, n = 0,..., N. 

Consider now the discretization scheme {Y^,V^, X^) discussed in Section]^ We use the approxi¬ 
mation (4.17) for the conditional expectations that have to be computed at each time step n. So, for 
every point iyi,Vn,k, Xn,j) xVn x (|4.17[) gives 




= E En.. iVn,k,Xn,j)Sn,k,jPh{Y O, b) k,j) 


(5.1) 


a,bG{d,u} 


where Sn,k,jPh denotes the operator in (4.18) applied to the function Ph{{n -|- l)/i, •), that is 


Sn,k,jPh{,^i b) 

Y’h ^(^ T Y)h,yi p iVn+lpaiPitk) ^n,A:) T P2y/'^^n,k iXn+l,jij(n,j) j ) j ^n-|-l,fca(n,fc) j • 

(5.2) 

We finally summarize the backward induction giving our approximating algorithm. For n = 0,1,..., A^, 
we define Ph{nh, y, v, x) for (y, v, x) G Ym x V)) x as follows: 

Ph{T, yi, VN,k, XN,j) = Y{yi) and as n = iV - 1,..., 0: 
pnh,yi,Vn,k,^nj) — 

- E En.< iYn,k j Xn,j )Pabi^^ 2 )'Sn,k,j (^) ^) ^ • 

a,b£{d,u} i&PM 

(5.3) 

Notice that, by (5.2), the computation of Sn,k,jPhi^,ci,b) requires the knowledge of the function 
y I—)• Ph{{n + l)/i, y, u,x) in points y’s that do not necessarily belong to the grid Ym- Therefore, in 
practice we compute such a function by means of quadratic interpolations. 
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Remark 5.1 Let us stress that the r.h.s. of (5.1) can be read in two equivalents ways. First, the 
term 

{Vn,k,Xn,j)Sn,k,jPhi^, a, 6), o, 6 £ {d, u}, i £ 3^m, 


is the numerical solution to the PDF (4.8) with final condition as in (5.2), so the r.h.s. of (|5.1|) is 


actually a weighted sum of the four solutions from each jump node {a, b) £ {d, u} for the pair {V, X), 
with weights given by the jump probabilities. But since the differential operator is linear in the Cauchy 
conditions, then one can first do the weighted sum of the final conditions, that is 

'Sn,k,jPh{^,a,b)p’fffn,k,j), i^yM, 

a,b£{d,u} 


and then apply the matrix Il{vn,k,Xn,j), i.e. solve the PDF (4.8) just once, and this is of course 
computationally less expensive. 

We can resume the main steps of our algorithm as follows. 

• Preprocessing: 


set the lattice Xn,j, 0 < j < n < iV, for the process X by using (3.1); 


set the lattice Vn^k-, 0 < fc < n < iV, for the process the V by using (3.5); 


merge the above lattices in a bivariate one {vn,k^Xn,j), 0 < k,j <n <N,hy using (3.9) 


— compute the jump-nodes and the transition probabilities Pab, (a, &) £ {d,u}, using (3.10); 

— set a mesh grid yt, i G yM, for the solution of all the PDE’s. 

Step N: for each node {vN,k,XN,j), 0 < k, j < N, compute the option prices at maturity for 
each Pi, i G by using the payoff function. 

Step n = N — 1,.. .0: for each {vn,k, Xn,j), 0 < k, j < n, compute the option prices for each pi, 
i G yM, by solving PDE (4.8) through ( |4.15 ), with terminal condition given by the weighted sum 
of the values at nodes (a, b) G {u, d} which have been computed in the previous step - weight by 


using the transition probabilities Pab (recall Remark 5.1). 


The theoretical proof of the convergence of our method is postponed to a further study. Although 
the ideas inspiring the method mainly come from [ 6 ] , here the convergence problem has to be tackled 
differently. In fact, in [ 6 ] the numerical scheme is written through a matrix 11 which is stochastic, so one 
can link the scheme to a Markov chain that approximates the process {Y, V, X) and use probabilistic 
methods (weak convergence) in order to study the convergence. But the scheme proposed here is 
purely numerical: the matrix n(u,x) = A~^{v,x) in (4.14) is stochastic if and only if /3 < |q:|, so the 


link with Markov chains fails and the probabilistic weak convergence cannot be used anymore. So, 
here we restrict ourselves to the study of the behavior and the efficiency of the proposed approach 
from the numerical point of view, see next Section 


6 Generalization to the Heston-Hull-White2d model 

The Heston-Hull-White2d model generalizes the previous model in the fact that the quantity g is 
assumed to be stochastic and to follow a diffusion model itself. So, the underlying process is now 
4-dimensional and is given by: the share price S, the volatility process V, the interest rate r and the 
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continuous dividend rate rj. Actually, here the process rj has not necessarily the meaning of a dividend 
rate, being for example a further interest rate process. In fact, the Heston-Hull-White2d model occurs 
in multi-currency models with short-rate interest rates, see e.g. m- 

Under the risk neutral measure, the dynamics are governed by the stochastic differential equation 

^ = in- r]t)dt + \/Vt dZt, 

dVt = Kv{0v - Vt)dt + avVVtdWl, 
drt = Kr{6r{t) — rt)dt + UrdW^, 
dr]t = - r]t)dt + ar^dWf, 

with initial data So,Vo,ro,'r]o > 0, where Z, W^, and denote possibly correlated Brownian 
motions. Note that the process r] evolves as a generalized OU process: Or^ is a deterministic function 
of the time. 

We consider non null correlations between the Brownian motions driving the pairs (S', U), (S, r) 
and [8,7]), that is 


d{Z, W^)t = Pi dt, d{Z, W^)t = P2 dt, d{Z, W^)t = P3 dt. 


Correlations among the processes V, r and p can be surely inserted (see next Remark 6.1). 

As done in Section]^ we take into account the transformations (2.1)-(2.2) for the generalized OU 
processes: we set 

rt = UrXl + (fl and pt = CTrjX^ + ip^ (6.1) 

where 


XI = -Kr / AJ ds + wf, pI = roe-^ 


/ er{s)e-'^^d-s)^g^ 

Jo Jo 

So ds + Bf, p^^ = poe-^^^ + K^ [ e^{s)e-^^d-^Us. 

Jo 


( 6 . 2 ) 


So, by considering the log-price process, we reduce to the 4-dimensional process (U, U, A^, A^) whose 
dynamics is given by 


dYt = priVt, XI, Xl t)dt + {pidWl + p2dWi + p^dWf + p^dW ^), 
dVt = pv{Vt)dt + avVVtdWS, 

dXl = px- {Xl)dt + dWS, (6.3) 

dX]' = pxv{XS)dt + dWf, 

with Yq = InSo G M, Uo > 0, A^ = 0, Aj = 0 


where 


with pI + pI + p^Y 1, 

Py{v, Xi, X2,t) = arXi -hPt - CrrtX2 ^ 

Pv{v) = Kv{dv — v), px^{x) = —KrX, px'7 (x) = — Kr/X. 

Starting from ( |6.3[ ) , we set-up an approximating procedure similar to the one developed in Section 
and Section In the following, we briefly describe how to extend such algorithms to the Heston- 
Hull-White2d model. 
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6.1 Approximation of (V, X’’, 

Concerning the triple we build an approximating tree on as follows: 


we apply the procedure in Section 3.1 to the process X^; 


we apply the procedure in Section 3.1 to the process X^-, 


• we apply the procedure in Section to the process V. 

We then get three approximating trees: 

Then, we use the null correlation between any two of V, X'’ and X^: we concatenate the above trees 
in order to get a 3-dimensional approximating tree (V^, X^’^, X^’^) for {V,X^,X'^) by introducing 
product-type jump probabilities. In other words, we generalize the probabilities in (3.10) for all the 
2^ = 8 possible jumps. 

Remark 6.1 One might include correlations between any two of the Brownian motions driving the 
processes V, X^ and X^. As described in Remark 3.1, the jump probabilities are no more of a 
product-type but they solve a linear system of equations that must include the matching of the local 
cross-moments up to order one in h. 


6.2 The scheme on the T-component and the approximating 4-dimensional process 


We repeat the reasonings in Section |4.1| in order to define an approximating time-continuous process 
{Y^,V^, X^'^, X^’^) for (Y,V,X^,X^) - roughly speaking, it suffices to replace the one-dimensional 
process X in Section 4.1 with the 2-dimensional process (X’’,X’?). So, we start from 

dYt = ^dVt + P2VVtdXl + psy^tdXf + p{Vt,Xf,Xf, t)dt + p^^t dWf (6.4) 

cry 

with 

p{v,Xi,X 2 ,t) = pY{v,Xi,X2,t) - —pv{v) - P 2 ^/vpx-{xi) - PsVv PX^ix2)- (6.5) 

(Tv 

Then, we apply the finite-difference method in Section |4.2| and we obtain a final difference scheme 
given by 

n(u, Xl,X2) = A~^{V, Xl,X2) 


where, p{-) being defined in (6.5) and A is given in (4.12) with 

a =-^ p{v,xi,X 2 ,nh) and ^ pIv■ 


( 6 . 6 ) 


Finally, we extend the approximation scheme (4.17) to the case in which X = (X’’,X’?) and the 
algorithm for the pricing of European or American options described in Section 


Remark 6.2 Let us briefly discuss the complexity of our algorithms. At each time step n = 0,..., N = 
T/h one has to find the solution of a PDE on a grid with 2M Y 1 points for eaeh fixed values of 

• case 1, Heston-Hull-White model: the pair (v,x) G x 

• case 2, Heston-Hull-White2d model: the triple {v,xi,X 2 ) G x x 


14 












The cardinality of all these possible values in case i is at most n x n^, i = 1,2. For each case, 
the system of equations (4.11) with tridiagonal matrix ( 4.12| ), can be solved by an efficient form of 
Gaussian elimination requiring a linear cost of order 0{M). Therefore, the total cost of our approach 
is of order 


N 

E 

n=l 


n 


i+1 


X (2M + 1) = 0(iV*+2 X M), case i = 1,2. 


We notice that the use of a full finite-difference scheme could be more expensive for practical compu¬ 
tations. Indeed, consider case 1 (Heston-Hull-White model). The solution of a 3-dimensional problem 
by applying finite-differences in all three components leads to the inversion of a big band-matrix. To 
reduce the computational cost, the problem requires to apply appropriate techniques such as ADI (Alter¬ 
nating Direction Implicit) techniques, see fW\/ and references therein. Specifically, in the authors 
propose an ADI approach to solve the Heston-Hull-White partial differential equation which needs a 
non-trivial implementation effort with a computational cost at least of order 0{M^) per time step, 
so the total cost is of order 0[N x M^) . Furthermore, as the dimension of the problem increases, 
it is not clear what happens if the problem is solved with a full finite-difference scheme. In case 2 
(Heston-Hull-White2d model), one should solve a 4-dimensional problem, bringing to the inversion of 
a very big band-matrix. This would give a cost which is hard to be quantified, and possibly in such a 
case the costs of the two procedures are no longer comparable. 


7 The hybrid Monte Carlo algorithm 


The approximation we have set-up for the Heston-Hull-White processes can be used to construct a 
Monte Carlo algorithm. Let us see how one can simulate a single path by using the tree approximation 
and the standard Euler scheme for the T-component. We call it “hybrid” because two different noise 
sources are considered: we simulate a continuous process in space (the component Y) starting from a 
discrete process in space (the 3-dimensional tree for {V, X'', X^)). 

Concerning the Heston-Hull-White dynamics in Section]^ consider the triple {Y,V,X) as in (2.3). 
Let {Vn , Xlf )n=o,i,...,N denote the Markov chain that approximates the pair {V,X). We construct 
a sequence {Yn)n=o,i,...,N approximating Y at times n = 0,1,..., by means of the Euler scheme 


defined in (4.3): we set Yff = Lq and for t G [nh, (n -|- l)h] with n = 0,1,..., — 1 then 


= yO + ^(^"+1 - Vjf) + P 2 V^(l ^+1 - Xd) + fiivO,Xf, nh)h + An+ 1 , (7.1) 

Uv 


where p is defined in ( |4.2[ ) and Ai,..., A^v denote i.i.d. standard normal r.v.’s, independent of the 
noise driving the chain {V,X). So, the simulation algorithm is very simple: at each time step n > 1, 


one let the pair {V,X) evolve on the tree and simulate the process Y at time nh by using (7.1). 


A similar algorithm can be considered to simulate the Heston-Hull-White2d dynamics in Section 
that can be seen as a function of the triple {Y,V, X^, X'^) in (6.3). Here, we apply the Euler 
scheme to ( [^ . So, let {Vj), Xf)^, X(l’’^)n=o,i,...,N denote the Markov chain approximating {V, X"^, A^), 


as described in Section 6.1, Starting from (6.4), we approximate the component Y at times nh, 
n = 0,1,..., A, as follows: we set Yq = Yq and for n = 1,..., A, n = 0,1,..., A — 1 then 




av 


'Vh(X^T _ vrM 




+P( 14 ^ nh)h + 


(7.2) 
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where ^ is defined in ( |6.5[ ) and Ai,..., A^v denote i.i.d. standard normal r.v.’s, independent of the 
noise driving the chain (V^, X^’^, And again, the simulation algorithm is straightforward. 


8 Numerical results 


In this section we provide numerical results in order to asses the efficiency and the robustness of our 
hybrid numerical approach. We first consider test experiments for the Heston-Hull-White model for 
the computation of European, American and barrier options (Section 8.1) and, following Andersen 
[3], we study Vanilla options with large maturities when the Feller condition is not fulfilled (Section 
|8.2|). Then we test European and American options in the Heston-Hull-White2d model (Section |8.3[). 


8.1 European, American and barrier options in the Heston-Hull-White model 

In the European and American option contracts we are dealing with, we consider the following set of 
parameters: 

• initial share price Sq = 100, strike price K = 100, maturity T = 1, dividend rate rj = 0.03; 

• initial interest rate tq = 0.04, speed of mean-reversion Kr = 1, interest rate volatility = 0.2, 
time-varying long-term mean 9r{t) which fits the theoretical bond prices to the yield curve 
observed on the market - to this purpose, we have chosen the interest rate curve given by 
P,(0,r) = e-0-0^^; 

• initial volatility Vq = 0 .1, long-mean Oy = 0.1, speed of mean-reversion Ky = 2, volatility of 
volatility uy = 0.3; 

• varying correlations: for the pairs {S,V), and {S,r), we set pi = psy = —0.5 and p 2 = psr = 
—0.5,0, 0.5 respectively; no correlation is assumed to exist between r and V. 


We notice that, under the above requests, the Feller condition holds. We postpone to next Section 


8.2 the analysis of cases in which the Feller condition is not fuffilled. 

The numerical study of the hybrid tree/finite-difference method HTFD is split in two cases: 


- HTFDl refers to the (fixed) number of time steps Nt = 50 and varying number of space steps 
Ns = 50,100,150,200; 


- HTFD2 refers to Nt = Ns = 50,100,150, 200. 


Concerning the Monte Carlo method, we compare the results by using the hybrid simulation scheme 
in Section that we call HMC. We also simulate paths by using the accurate third-order Alfonsi 
[I] discretization scheme for the CIR stochastic volatility process and by using an exact scheme for 
the interest rate. These simulating schemes are here called AMC. In both Monte Carlo methods, 
we consider varying number of time discretization steps Nt = 50,100,150, 200 and two cases for the 
number of Monte Carlo iterations: 


- HMCl and AMCl refer to 50 000 iterations, 

- HMC2 and AMC2 refer to 200 000 iterations. 
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In the European case, the benchmark value B-AMC is computed using the Alfonsi method with 
300 discretization time steps and the associated Monte Carlo estimator is computed with 1 million 
simulations. In the American case, in absence of reliable numerical methods, the benchmark values 
B-AMC-LS are obtained by the Longstaff-Schwartz m Monte Carlo algorithm with 50 exercise 
dates, combined with the Alfonsi method with 300 discretization time steps and 1 million iterations. 

Table reports both European call option prices and implied volatilities results. In Table we 
provide American call option prices. Table refers to the computational time cost (in seconds) of the 
different algorithms in the call European case. 

The numerical results show that HTFD is accurate, reliable and efficient for pricing European 
and American options in the Heston-Hull-White model. Moreover, our hybrid Monte Carlo algorithm 
HMC appears to be competitive with AMC, that is the one from the accurate simulations by Alfonsi 
[1]: the numerical results are similar in term of precision and variance but HMC is definitely better 
from the computational times point of view. Additionally, because of its simplicity, HMC represents 
a real and interesting alternative to AMC. As a further evidence of the accuracy of our methods, in 
Eigure 2 we study the shapes of implied volatility smiles across moneyness ^ using HTFDl with 
Nt = 5(J and Ns = 200 and HMCl with Nt = 50, and we compare the graphs with the results from 
the benchmark. 

In order to study the convergence behavior of our approach HTFD, we consider the convergence 
ratio proposed in m, defined as 

Pn — Pn 

I'atio = (8.1) 

2 

where P^ denotes here the approximated price obtained with N = Nt number of time steps. Recall 
that Pn = 0{N~°^) means that ratio = 2“. For the sake of comparison with the numerical convergence 
speed studied in [6], we report ratios for American put options. We split the numerical study in two 
different cases: when the Feller condition holds and when it does not, the results being given in 
Table 1^ and Table respectively (details on the option parameters are given in the table captions). 
Both tables give evidence of the numerical convergence, but with some differences. In fact, under the 
Feller condition (Table Q , the numerical speed of convergence is definitely linear (this is not really 
surprising because tree methods are usually linear), whereas in the opposite case (Tablethe behavior 
is approximately linear. 

Furthermore, we study the behavior of HTFD in the case of exotic options, namely for continu¬ 
ously monitored barrier options. We consider call up-and-out options, whose payoff is given by 

{St - K)Tl^St<Hvt<T}- 

In our numerical experiments, the up barrier is set at = 130 and we choose different values for 
S'o = 80,100,120. Table 1^ reports European call up-and-out option prices. In the barrier option case, 
we compare with a benchmark value, called B-AMC, computed by 2 millions iterations which use 
the Alfonsi AMC method with 9600 discretization time steps. The numerical results confirm the 
reliability of HTFD for barrier options. 
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(a) 




Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

PSr = 

-0.5 

50 

11.202744 

11.202744 

11.34±0.04 

11.30±0.16 

11.32±0.08 

11.34±0.16 

11.37±0.08 



100 

11.319814 

11.331040 


11.41±0.16 

11.38±0.08 

11.31±0.16 

11.36±0.08 



150 

11.340665 

11.349902 


11.36±0.16 

11.36±0.08 

11.35±0.16 

11.38±0.08 



200 

11.346972 

11.355772 


11.34±0.16 

11.37±0.08 

11.44±0.16 

11.39±0.08 

PSr 

= 0 

50 

12.526779 

12.526779 

12.77±0.04 

12.66±0.18 

12.69±0.09 

12.68±0.18 

12.79±0.09 



100 

12.720651 

12.705772 


12.74±0.18 

12.79±0.09 

12.63±0.18 

12.78±0.09 



150 

12.754610 

12.749526 


12.74±0.18 

12.79±0.09 

12.68±0.18 

12.81±0.09 



200 

12.760365 

12.766836 


12.74±0.18 

12.80±0.09 

12.75±0.18 

12.79±0.09 

PSr = 

= 0.5 

50 

13.853193 

13.853193 

14.04±0.04 

13.88±0.19 

13.92±0.10 

13.97±0.20 

14.05±0.10 



100 

14.011537 

14.013063 


13.91±0.19 

14.01±0.10 

13.89±0.19 

14.06±0.10 



150 

14.031598 

14.038361 


13.94±0.19 

14.07±0.10 

13.92±0.20 

14.08±0.10 



200 

14.038235 

14.045612 


13.99±0.19 

14.07±0.10 

13.90±0.19 

14.06±0.10 


(b) 




Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

PSr = 

-0.5 

50 

0.279002 

0.279002 

0.282602 

0.281649 

0.282117 

0.282602 

0.283389 



100 

0.282073 

0.282367 


0.284443 

0.283681 

0.281815 

0.283127 



150 

0.282620 

0.282862 


0.283034 

0.283085 

0.282865 

0.283652 



200 

0.282785 

0.283016 


0.282478 

0.283408 

0.285226 

0.283914 

PSr 

= 0 

50 

0.313772 

0.313772 

0.320169 

0.317398 

0.317958 

0.317802 

0.320695 



100 

0.318871 

0.318480 


0.319306 

0.320650 

0.316487 

0.320432 



150 

0.319764 

0.319630 


0.319063 

0.320716 

0.317802 

0.321221 



200 

0.319916 

0.320086 


0.319288 

0.321009 

0.319643 

0.320695 

PSr = 

= 0.5 

50 

0.348697 

0.348697 

0.353623 

0.349329 

0.350359 

0.351777 

0.353887 



100 

0.352873 

0.352913 


0.350234 

0.352954 

0.349667 

0.354151 



150 

0.353402 

0.353580 


0.350960 

0.354324 

0.350458 

0.354679 



200 

0.353577 

0.353771 


0.352184 

0.354545 

0.349931 

0.354151 


Table 1: Prices (a) and Implied volatilities (b) of European call options. Sq = 100, K = 100, T = 1, vq = 0.04, 
Kr = 1, ar = 0.2, rj = 0.03, Vq = 0-1, By = 0.1, Kv = 2, ay = 0.3, psr = —0.5, 0, 0.5, psy = —0.5. 



Ns 

HTFDl 

HTFD2 

B-AMC-LS 

pSr = -0.5 

50 

12.090433 

12.090433 

12.22±0.01 


100 

12.205014 

12.212884 



150 

12.224432 

12.231392 



200 

12.230288 

12.237054 


pSr = 0 

50 

12.912708 

12.912708 

13.16±0.02 


100 

13.119121 

13.101073 



150 

13.156492 

13.149182 



200 

13.162893 

13.168602 


pSr ~ 0.5 

50 

13.944266 

13.944266 

14.15±0.02 


100 

14.125059 

14.122918 



150 

14.146240 

14.152060 



200 

14.153288 

14.160288 



Table 2: Prices of American call options. Sq = 100, K = 100, T = 1, ro = 0.04, = 1; = 0.2, rj = 0.03, 

Vo = 0.1, By = 0.1, Ky = 2, ay = 0.3, psr = —0.5, 0, 0.5, Psy = —0.5. 
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Ns 

HTFDl 

HTDF2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC 

50 

0.41 

0.41 

223.67 

0.77 

3.05 

2.16 

7.48 

100 

0.84 

11.33 


1.59 

6.11 

4.00 

14.61 

150 

1.37 

49.99 


2.33 

9.13 

5.87 

21.64 

200 

1.87 

213.06 


3.11 

12.73 

7.61 

28.85 


Table 3: Computational times (in seconds) for European call options. 


Inplied Volatility Snile 



Figure 2: Moneyness vs implied volatility for European call options. T = 1, tq = 0.04, Kr = 1, 
Ur = 0.2, rj = 0.03, Vo = 0-1, Oy = 0.1, Ky = 2, ay = 0.3, psr = —0.5, psy = —0.5. 


K 

Nt 

N 

Price 

Ratio 

80 

25 

50 

21.494606 



50 

100 

21.534555 



100 

200 

21.553473 

2.111762 


200 

400 

21.563911 

1.812303 


400 

800 

21.569080 

2.019428 

100 

25 

50 

12.607035 



50 

100 

12.749006 



100 

200 

12.815657 

2.130053 


200 

400 

12.845050 

2.267634 


400 

800 

12.859561 

2.025489 

120 

25 

50 

21.444819 



50 

100 

21.539534 



100 

200 

21.572106 

2.907912 


200 

400 

21.586338 

2.288708 


400 

800 

21.592706 

2.234825 


Table 4; HTFD-ratio (8.1) for the price of American put options at final time T = 0.25. Sq = 100, p = 0.03, 
To = 0.04,fcr = 1, CTr = 0.2 , Vq = 0 -1, ky = 2, Oy = 0.1, (Ty = 0.3, psv = —0.5, pSr = 0.5. 
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K 

Nt 

N 

Price 

Ratio 

80 

25 

50 

21.635830 



50 

100 

21.669504 



100 

200 

21.688879 

1.738049 


200 

400 

21.700965 

1.603169 


400 

800 

21.710373 

1.284610 

100 

25 

50 

10.649104 



50 

100 

10.762867 



100 

200 

10.812709 

2.282480 


200 

400 

10.835512 

2.185787 


400 

800 

10.848349 

1.776369 

120 

25 

50 

20.755654 



50 

100 

20.873859 



100 

200 

20.908825 

3.380584 


200 

400 

20.919694 

3.216994 


400 

800 

20.924295 

2.362300 


Table 5: HTFD-ratio 

To = 0.0A,kr = 1, ar = 


( 8 . 1 ) for the price 
0.2, To = 0.09, kv 


of American put options at final time T = 0.25. Sq = 100, 77 = 0.03, 
= 1, Ov — 0.09, ay = 1, psv = —0.3, psr = 0. 




Ns 

HTFDl 

HTFD2 

B-AMC 

So = 

: 80 

50 

1.211544 

1.211544 




100 

1.251453 

1.255849 

1.282211±0.01 



150 

1.264327 

1.270193 




200 

1.269703 

1.274332 


So = 

100 

50 

1.819848 

1.819848 




100 

1.941320 

1.916440 

1.947565±0.01 



150 

1.964666 

1.930681 




200 

1.974201 

1.933482 


So = 

120 

50 

0.697718 

0.697718 




100 

0.749116 

0.725243 

0.728431±0.01 



150 

0.762224 

0.726872 




200 

0.766022 

0.725139 



Table 6: Prices of European call up-and-out options. Up harrier is H = 130. K = 100, , T = 1, tq = 0.04, 
Kr = 1, CTr = 0.2, T] = 0.03, Vq = 0.1, Oy = 0.1, Ky = 2, ay = 0.3, psr = —0.5, psy = —0.5. 
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8.2 European options with large maturity in the Heston-Hull-White model 

In order to verify the robustness of the proposed algorithms we consider experiments when the Feller 
condition is not fulfilled and with large maturities. We test here the cases I, II, III (reordered with 
respect to the maturity) proposed in Andersen [H] in order to price European call options. Moreover, 
we add the case IV with maturity T = 25. 

We consider the following values for the parameters of the model and for the maturity date: 


Case I: Vq = 0.09, 9y = 0.09, Ky = 1, o'v — 1? Psv — ~0.3, T = 5; 

Case II: Vq = 0.04, 9y = 0.04, Ky = 0.5, ay = 1, psy = —0.9, T = 10. 
Case III: Vq = 0.04, 9y = 0.04, Ky = 0.3, ay = 0.9, psy = -0.5, T = 15. 
Case IV: Vq = 0.04, 9y = 0.04, Ky = 0.3, ay = 0.9, psy = -0.5, T = 25. 


We take into account varying strikes K = 70,100,140. No correlation is assumed to exist between 
S and r, that is psr = 0, so we can compare the results with the semi closed-form analytic formula 
(SCF) for European call options which is available in |13j . We use in particular the implementation 
of the semi closed-form analytic formula provided in QuantLib |21j . Moreover in all cases the interest 


rate parameters, the initial share value and the dividend are the same of Section 8.1 
• So = 100, p = 0.03; 


• ro = 0.04, Kr = 1, ar = 0.2. 

In Tables 0 § [9, [To| we provide European call option prices and implied volatility results. The 
numerical results suggest that large maturities bring to a slight loss of accuracy for both HTFD and 
HMC, even if each method provides a satisfactory approximation of the true option prices. It is worth 
noticing that for long maturities T = 5,15, 25 we have developed experiments with the same number 
of steps both in time (W) and space (Ns) as for T = 1. So, the numerical experiments are not slower, 
and it is clear that one could achieve a better accuracy for larger values of Nt- 
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(a) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

37.054163 

37.054163 

37.491811 

37.36±0.47 

37.32±0.23 

37.38±0.47 

37.31±0.23 



100 

37.392491 

37.395372 


37.30±0.45 

37.52±0.24 

37.61±0.46 

37.61±0.24 



150 

37.480467 

37.521733 


37.40±0.46 

37.58±0.24 

37.55±0.47 

37.58±0.24 



200 

37.546885 

37.570675 


37.42±0.46 

37.48±0.23 

37.49±0.51 

37.60±0.24 

K = 

100 

50 

23.997806 

23.997806 

24.706195 

24.64±0.43 

24.58±0.21 

24.61±0.43 

24.54±0.21 



100 

24.537750 

24.540987 


24.49±0.41 

24.76±0.21 

24.79±0.43 

24.81±0.22 



150 

24.669356 

24.684708 


24.60±0.41 

24.81±0.22 

24.71±0.42 

24.78±0.22 



200 

24.747161 

24.766840 


24.67±0.42 

24.70±0.21 

24.73±0.47 

24.82±0.22 

K = 

140 

50 

13.672435 

13.672435 

14.324566 

14.33±0.38 

14.24±0.18 

14.22±0.37 

14.17±0.18 



100 

14.248533 

14.205762 


14.11±0.35 

14.40±0.198 

14.40±0.37 

14.40±0.19 



150 

14.373163 

14.318446 


14.21±0.36 

14.43±0.198 

14.33±0.38 

14.40±0.19 



200 

14.444183 

14.404071 


14.31±0.36 

14.32±0.188 

14.31±0.42 

14.42±0.20 


(b) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

0.313372 

0.313372 

0.322137 

0.319432 

0.318614 

0.319863 

0.318541 



100 

0.320152 

0.320209 


0.318384 

0.322782 

0.324521 

0.324567 



150 

0.321910 

0.322734 


0.320315 

0.323861 

0.323277 

0.323956 



200 

0.323236 

0.323711 


0.320737 

0.321815 

0.322027 

0.324318 

K = 

100 

50 

0.296912 

0.296912 

0.306954 

0.306002 

0.305124 

0.30564 

0.304539 



100 

0.304563 

0.304608 


0.303947 

0.307727 

0.308148 

0.308367 



150 

0.306431 

0.306649 


0.305385 

0.308440 

0.307022 

0.307959 



200 

0.307536 

0.307815 


0.306431 

0.306889 

0.307262 

0.308556 

K = 

140 

50 

0.291198 

0.291198 

0.299737 

0.299844 

0.298690 

0.298395 

0.240057 



100 

0.298743 

0.298183 


0.296939 

0.300702 

0.240301 

0.239857 



150 

0.300373 

0.299657 


0.298282 

0.301138 

0.299848 

0.300773 



200 

0.301301 

0.300777 


0.299533 

0.299736 

0.299505 

0.301033 


Table 7; Prices (a) and Implied volatilities (b) of European call options. So = 100, T = 5, tq = 0.04, Kr = 1, 
(Tr = 0.2, rj = 0.03, Vo = 0.09, 0v = 0.09, kv = 1; cry = 1, psr = 0, psv = —0.3, K = 70,100,140. 
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(a) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

33.702753 

33.702753 

34.101622 

33.59i0.23 

33.76i0.ll 

34.12i0.23 

34.09i0.ll 



100 

33.773407 

34.120510 


33.98i0.23 

34.10i0.ll 

34.25i0.23 

34.10i0.ll 



150 

33.776196 

33.818752 


33.61i0.23 

33.76i0.ll 

34.02i0.23 

34.lli0.ll 



200 

33.778268 

33.944743 


33.91i0.23 

33.91i0.ll 

34.00i0.23 

34.10i0.ll 

K = 

100 

50 

22.540546 

22.540546 

23.140518 

22.57i0.21 

22.65i0.10 

23.14i0.21 

23.09i0.ll 



100 

22.761622 

23.076646 


22.95i0.21 

23.06i0.10 

23.26i0.21 

23.12i0.ll 



150 

22.795766 

22.857113 


22.68i0.21 

22.81i0.10 

23.04i0.21 

23.09i0.ll 



200 

22.806087 

22.978809 


22.96i0.21 

22.95i0.10 

23.02i0.21 

23.13i0.ll 

K = 

140 

50 

13.335726 

13.335726 

13.755466 

13.18i0.17 

13.21i0.08 

13.72i0.17 

13.68i0.09 



100 

13.510432 

13.726749 


13.53i0.17 

13.62i0.09 

13.86i0.18 

13.72i0.09 



150 

13.528322 

13.553294 


13.34i0.17 

13.46i0.09 

13.66i0.17 

13.73i0.09 



200 

13.530288 

13.639595 


13.60i0.17 

13.60i0.09 

13.64i0.17 

13.76i0.09 


(b) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

0.227844 

0.227844 

0.234811 

0.225850 

0.228795 

0.235187 

0.234676 



100 

0.229082 

0.235140 


0.232714 

0.234866 

0.237332 

0.234768 



150 

0.229131 

0.229876 


0.226245 

0.228809 

0.233392 

0.235042 



200 

0.229167 

0.232077 


0.231443 

0.231474 

0.232965 

0.234723 

K = 

100 

50 

0.215548 

0.215548 

0.222789 

0.215951 

0.216908 

0.222801 

0.222156 



100 

0.218214 

0.222017 


0.220435 

0.221799 

0.224184 

0.222572 



150 

0.218625 

0.219366 


0.217216 

0.218739 

0.221545 

0.222656 



200 

0.218750 

0.220835 


0.220583 

0.220512 

0.221377 

0.222760 

K = 

140 

50 

0.210662 

0.210662 

0.215154 

0.209030 

0.209283 

0.214777 

0.214386 



100 

0.212532 

0.214847 


0.212764 

0.213679 

0.216253 

0.214726 



150 

0.212723 

0.212991 


0.210669 

0.212018 

0.214162 

0.214846 



200 

0.212744 

0.213914 


0.213460 

0.213506 

0.213908 

0.215215 


Table 8 ; Prices (a) and Implied volatilities (b) of European call options. So = 100, T = 10, ro = 0.04, Kr = 1, 
cr^ = 0.2, p = 0.03, Vo = 0.04, 0v = 0.04, kv = 0.5, cry = 1, psr = 0, psv = -0.9, K = 70,100,140. 
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(a) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

32.872766 

32.872766 

33.182814 

33.17±0.31 

33.26±0.16 

33.13dz0.31 

33.18zb0.16 



100 

33.041266 

33.161213 


33.10±0.30 

33.29±0.15 

33.18dz0.31 

33.19zb0.15 



150 

33.098186 

33.159078 


33.03±0.30 

33.12±0.16 

33.20zb0.34 

33.25zb0.16 



200 

33.150052 

33.235555 


33.02±0.29 

33.11±0.15 

33.12zb0.33 

33.35zb0.15 

K = 

100 

50 

24.738008 

24.738008 

25.183109 

25.00±0.30 

25.05±0.15 

25.10zb0.30 

25.17zb0.15 



100 

24.979024 

25.089961 


24.96±0.29 

25.18±0.15 

25.20zb0.30 

25.20zb0.15 



150 

25.047214 

25.150207 


24.99±0.28 

25.11±0.15 

25.17zb0.33 

25.23zb0.16 



200 

25.103492 

25.224136 


24.97±0.28 

25.09±0.15 

25.07zb0.31 

25.30zb0.15 

K = 

140 

50 

17.522401 

17.522401 

17.851374 

17.49±0.27 

17.53±0.14 

17.76zb0.27 

17.84zb0.14 



100 

17.702990 

17.779408 


17.51±0.26 

17.74±0.14 

17.85zb0.28 

17.86zb0.15 



150 

17.752550 

17.858103 


17.59±0.26 

17.77±0.14 

17.82±0.31 

17.87zb0.14 



200 

17.800293 

17.912261 


17.59±0.25 

17.75±0.13 

17.73zb0.29 

17.93zb0.14 


(b) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

0.231761 

0.231761 

0.237013 

0.236812 

0.238369 

0.236053 

0.236928 



100 

0.234617 

0.236648 


0.235577 

0.238826 

0.236974 

0.237216 



150 

0.235581 

0.236612 


0.234478 

0.235946 

0.237321 

0.238081 



200 

0.236459 

0.237905 


0.234214 

0.235730 

0.235877 

0.239816 

K = 

100 

50 

0.225390 

0.225390 

0.230837 

0.228633 

0.22926 

0.229786 

0.230708 



100 

0.228336 

0.229695 


0.228045 

0.230802 

0.231012 

0.231068 



150 

0.229171 

0.230433 


0.228424 

0.229935 

0.230683 

0.231375 



200 

0.229861 

0.231340 


0.228236 

0.229682 

0.229474 

0.232223 

K = 

140 

50 

0.223601 

0.223601 

0.227024 

0.223225 

0.223688 

0.226081 

0.226894 



100 

0.225479 

0.226275 


0.223424 

0.225813 

0.227054 

0.227136 



150 

0.225995 

0.227094 


0.224321 

0.226215 

0.226728 

0.227195 



200 

0.226492 

0.227659 


0.224316 

0.225961 

0.225781 

0.227881 


Table 9; Prices (a) and Implied volatilities (b) of European call options. So = 100, T = 15, ro = 0.04, Kr = 1, 
cr^ = 0.2, p = 0.03, Vo = 0.04, 0v = 0.04, kv = 0.3, cry = 0.9, psr = 0, psv = -0.5, K = 70,100,140. 
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(a) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

: 70 

50 

28.772135 

28.772135 

28.969593 

29.01±0.29 

29.05±0.15 

28.92±0.29 

28.98±0.15 



100 

28.890859 

29.076376 


29.06±0.34 

29.00±0.15 

28.99±0.29 

28.97±0.15 



150 

29.007171 

29.225059 


29.07±0.29 

29.15±0.15 

28.90±0.31 

29.05±0.15 



200 

29.125812 

29.152251 


29.05±0.31 

28.91±0.15 

28.95±0.29 

29.01±0.14 

K = 

100 

50 

23.947048 

23.947048 

24.255944 

24.09±0.28 

24.13±0.15 

24.20±0.29 

24.26±0.15 



100 

24.107443 

24.300298 


24.22±0.33 

24.19±0.155 

24.27±0.28 

24.25±0.15 



150 

24.233382 

24.462163 


24.26±0.28 

24.37±0.155 

24.17±0.31 

24.33±0.15 



200 

24.356051 

24.436578 


24.32±0.31 

24.20±0.145 

24.22±0.28 

24.31±0.14 

K = 

140 

50 

19.352114 

19.352114 

19.601699 

19.21±0.27 

19.24±0.14 

19.52±0.27 

19.59±0.14 



100 

19.459177 

19.637550 


19.39±0.32 

19.39±0.14 

19.62±0.27 

19.59±0.14 



150 

19.567765 

19.778396 


19.51±0.27 

19.62±0.14 

19.51±0.30 

19.66±0.14 



200 

19.692584 

19.798050 


19.60±0.30 

19.53±0.14 

19.55±0.27 

19.65±0.13 


(b) 




Ns 

HTFDl 

HTFD2 

SCF 

HMCl 

HMC2 

AMCl 

AMC2 

K = 

^ 70 

50 

0.235972 

0.235972 

0.239830 

0.240620 

0.241466 

0.238797 

0.240057 



100 

0.238291 

0.241919 


0.241551 

0.240401 

0.240301 

0.239857 



150 

0.240565 

0.244831 


0.241698 

0.243433 

0.239857 

0.241365 



200 

0.242887 

0.243405 


0.241403 

0.238618 

0.239442 

0.239442 

K = 

100 

50 

0.231127 

0.231127 

0.235633 

0.233270 

0.233862 

0.234826 

0.235633 



100 

0.233464 

0.236283 


0.235133 

0.234655 

0.235771 

0.235536 



150 

0.235303 

0.238656 


0.235636 

0.237270 

0.234367 

0.236650 



200 

0.237099 

0.238281 


0.236537 

0.234765 

0.235096 

0.236425 

K = 

140 

50 

0.229635 

0.229635 

0.232641 

0.227954 

0.228325 

0.231645 

0.232511 



100 

0.230923 

0.233074 


0.230084 

0.230128 

0.233016 

0.232455 



150 

0.2332231 

0.234777 


0.231504 

0.232913 

0.231583 

0.233359 



200 

0.2333739 

0.235015 


0.232634 

0.231717 

0.232046 

0.233281 


Table 10: Prices (a) and Implied volatilities (b) of European eall options. So = 100, T = 25, rg = 0.04, Kr = 1, 
cr^ = 0.2, p = 0.03, Vo = 0.04, 0v = 0.04, kv = 0.3, cry = 0.9, psr = 0, psv = -0.5, K = 70,100,140. 
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8.3 European and American options in the Heston-Hull-White2d model 

In the European and American option contracts we are dealing with, we consider the following set of 
parameters: 


• So = 100, K = 100, T = 1; 

• ro = 0.04, r]o = 0.03, Kr = Krj = 1, ar = = 0.2; 

• Vo = 0.1, 9y = 0.1, Ky = 2, (Ty = 0.3; 

• PSr — 0.5,0, 0.5, PSV — 0.5, PSri — 0.5, 0.5, pyr — PVr) — Prrj — 0, 

• Pr{Q,T) = P^(0,r) = e-o-03'r. 


As before, the time-varying long-term means 6r{t) and 6n{t) fit the theoretical bond prices Pr{0,T) 
and P,j(0,T) to the yield curve observed on the market. We make this choice following the multi- 
currency models with short-rate interest rates in [16]. We consider here only the number of space 
steps Ns = 30, 50,100 because the cases Ns = 150, 200 need a too high computational time. Tables 
11 [12 and |13| report European and American call option prices and implied volatilities. As before. 


the benchmark value for European options is computed using the Alfonsi B-AMC method with 
300 discretization time steps and the associated Monte Carlo estimator is computed with 1 million 
iterations. Concerning the benchmark B-AMC-LS for American options, it is computed by means 
of the Longstaff-Schwartz m Monte Carlo algorithm with 50 exercise dates, combined with the 
Alfonsi method with 300 discretization time steps and 1 million iterations. Table 14 refers to the 
computational time cost (in seconds) of the different algorithms in the call European case. In Figure]^ 
we compare the shapes of implied volatility smiles across moneyness ^ using HTFDl with Nt = 30 
and Ns = 100 and HMCl with Nt = 30. The numerical results confirm the good numerical behavior 
of HTFD and HMC in the Heston-Hull-White2d model as well. 
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(a) 


psv = 

PSr, = 

-0.5, 

-0.5 

Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

PSr = 

-0.5 

30 

13.470572 

13.470572 

13.79 ± 0.04 

13.82±0.20 

13.74±0.10 

13.83±0.20 

13.79±0.10 



50 

13.688842 

13.671173 


13.96±0.20 

13.81±0.10 

13.88±0.20 

13.80±0.10 



100 

13.790205 

13.781519 


14.00±0.20 

13.80±0.10 

13.68±0.20 

13.73±0.10 

PSr 

= 0 

30 

14.736242 

14.736242 

15.04 ± 0.05 

15.10±0.22 

14.99±0.11 

14.95±0.22 

15.03±0.11 



50 

14.958094 

14.946029 


15.23±0.22 

15.04±0.11 

14.98±0.22 

15.01±0.11 



100 

15.019204 

15.032709 


15.21±0.22 

15.04±0.11 

14.80±0.21 

14.97±0.11 

PSr = 

= 0.5 

30 

15.805046 

15.805046 

16.19 ± 0.03 

16.13±0.23 

16.06±0.11 

16.04±0.23 

16.17±0.12 



50 

16.052315 

16.032043 


16.33±0.23 

16.10±0.11 

16.09±0.23 

16.13±0.12 



100 

16.155354 

16.145308 


16.24±0.23 

16.19±0.12 

15.93±0.23 

16.12±0.12 


(b) 


psv = —0.5, 
pSr, = -0.5 

Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

pSr = —0.5 

30 

0.338612 

0.338612 

0.347031 

0.347724 

0.345593 

0.348085 

0.347031 


50 

0.344364 

0.343898 


0.351511 

0.347675 

0.349404 

0.347294 


100 

0.347036 

0.346807 


0.352510 

0.347205 

0.344131 

0.345449 

pSr = 0 

30 

0.372004 

0.372004 

0.380033 

0.381610 

0.378689 

0.377653 

0.379769 


50 

0.377867 

0.377549 


0.385153 

0.380032 

0.378447 

0.379240 


100 

0.379483 

0.379840 


0.384419 

0.380061 

0.373689 

0.378182 

pSr = 0.5 

30 

0.400281 

0.400281 

0.410485 

0.408889 

0.407043 

0.406508 

0.409954 


50 

0.406834 

0.406297 


0.414273 

0.408039 

0.407833 

0.408894 


100 

0.409566 

0.409300 


0.411792 

0.410506 

0.403592 

0.408629 


Table 11: Prices (a) and Implied volatilities (b) of European call options. So = 100, K = 100, T = 1, ro = 0.04, 
Kr = 1) CTj. = 0.2, ryo = 0.03, k,, = 1, arj = 0.2, Vq = 0.1, 0v = 0.1, Ky = 2, ay = 0.3, psr = —0.5,0, 0.5, 
Psv = -0.5, psr, = -0.5. 
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(a) 


psv = 

PSr, = 

-0.5, 

0.5 

Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

PSr = 

-0.5 

30 

9.418513 

9.418513 

9.61 ± 0.03 

9.57±0.13 

9.62±0.07 

9.64±0.13 

9.66±0.07 



50 

9.552565 

9.532194 


9.57±0.13 

9.61±0.07 

9.65±0.13 

9.66±0.07 



100 

9.633716 

9.607339 


9.66±0.13 

9.62±0.07 

9.63±0.13 

9.63±0.07 

PSr 

= 0 

30 

10.916753 

10.916753 

11.18 ± 0.03 

11.15±0.15 

11.16±0.08 

11.07±0.15 

11.22±0.08 



50 

11.117050 

11.100343 


11.18±0.15 

11.16±0.08 

11.14±0.15 

11.22±0.08 



100 

11.178119 

11.173631 


11.16±0.15 

11.18±0.08 

11.08±0.15 

11.20±0.08 

PSr - 

= 0.5 

30 

12.203271 

12.203271 

12.55 ± 0.04 

12.44±0.17 

12.43±0.09 

12.47±0.17 

12.60±0.09 



50 

12.443197 

12.411406 


12.54±0.17 

12.44±0.09 

12.53±0.17 

12.59±0.09 



100 

12.552842 

12.522237 


12.45±0.17 

12.55±0.09 

12.45±0.17 

12.58±0.09 


(b) 


psv = 

PSr, = 

-0.5, 

0.5 

Ns 

HTFDl 

HTFD2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

PSr = 

-0.5 

30 

0.232267 

0.232267 

0.237277 

0.236285 

0.237525 

0.238062 

0.238586 



50 

0.235774 

0.235241 


0.236157 

0.237263 

0.238324 

0.238586 



100 

0.237898 

0.237208 


0.238622 

0.237412 

0.237800 

0.237800 

PSr 

= 0 

30 

0.271502 

0.271502 

0.278405 

0.277704 

0.277855 

0.275520 

0.279454 



50 

0.276754 

0.276316 


0.278277 

0.277937 

0.277356 

0.279454 



100 

0.278356 

0.278238 


0.277841 

0.278435 

0.275782 

0.278930 

PSr ~ 

= 0.5 

30 

0.305269 

0.305269 

0.314383 

0.311364 

0.311103 

0.312279 

0.315698 



50 

0.311575 

0.310739 


0.313992 

0.311570 

0.313857 

0.315435 



100 

0.314458 

0.313653 


0.311665 

0.314463 

0.311754 

0.315172 


Table 12: Prices (a) and Implied volatilities (b) of European call options. So = 100, K = 100, T = 1, ro = 0.04, 
Kr = 1) CTj. = 0.2, r]o = 0.03, k,, = 1, arj = 0.2, Vq = 0.1, 0v = 0.1, Ky = 2, cry = 0.3, psr = —0.5,0, 0.5, 
PSV = -0.5, psr, = 0.5. 
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(a) 


Psv = 

PSr, = 

-0.5, 

-0.5 

Ns 

HTFDl 

HTFD2 

B-AMC-LS 

PSr = 

-0.5 

30 

14.057963 

14.057963 

14.40 ± 0.02 



50 

14.290597 

14.263254 




100 

14.400377 

14.381552 


PSr 

= 0 

30 

14.989844 

14.989844 

15.32 ± 0.02 



50 

15.253011 

15.229151 




100 

15.320569 

15.331744 


PSr - 

= 0.5 

30 

15.826696 

15.826696 

16.28 ± 0.02 



50 

16.146080 

16.111559 




100 

16.270439 

16.248656 


(b) 

psv = 

PSn = 

-0.5, 

0.5 

Ns 

HTFDl 

HTFD2 

B-AMC-LS 

PSr = 

-0.5 

30 

11.598655 

11.598655 

11.72 ± 0.02 



50 

11.707669 

11.681873 




100 

11.775632 

11.743388 


PSr 

= 0 

30 

12.400256 

12.400256 

12.60 ± 0.02 



50 

12.579124 

12.561214 




100 

12.634969 

12.629401 


PSr - 

= 0.5 

30 

13.137621 

13.137621 

13.47 ± 0.02 



50 

13.380571 

13.341882 




100 

13.497053 

13.459978 



Table 13: Prices of American call options. Sq = 100, K = 100, T = 1, tq = 0.04, Kr = 1, (Jr = 0.2, rjo = 0.03, 
Krj = 1, an = 0.2, Vo = 0.1, 9v = 0.1, Ky = 2, cry = 0.3, psr = —0.5,0,0.5, psv = —0.5, psn = —0.5,0.5. 


Ns 

HTFDl 

HTDF2 

B-AMC 

HMCl 

HMC2 

AMCl 

AMC2 

30 

2.22 

2.22 

284.84 

0.60 

2.61 

1.79 

6.03 

50 

4.15 

24.56 


1.14 

4.19 

2.73 

9.58 

100 

7.95 

998.1 


2.02 

8.06 

5.05 

18.70 


Table 14: Computational times (in seconds) for European call options. 
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Inplied Volatility Snile 



Figure 3: Moneyness vs implied volatility for European call options. T = 1, ro = 0.04, Kr = 1, 
ar = 0.2, rjo = 0.03, = 1, cr,, = 0.2, Vq = 0.1, 9v = 0.1, Ky = 2, ay = 0.3, psr = —0.5, psy = —0.5, 

PSr^ = -0.5. 


9 Conclusions 

We have introduced a new hybrid tree/finite-difFerence method and a new Monte Carlo method for 
numerically pricing options in a stochastic volatility framework with stochastic interest rates. The 
numerical comparisons show that our methods provide a good approximation of the option prices with 
efficient time computations. 

Acknowledgements. The authors wish to thank Andrea Molent for useful remarks and for having 
implemented the Alfonsi Monte Carlo scheme and the Longstaff-Schwarz algorithm. 


References 

[1] A. Alfonsi (2010): High order discretization schemes for the CIR process: application to affine term 
structure and Heston models, Mathematics of Computation, 79, 209-237. 

[2] K. Amin, A. Khanna (1994): Convergence of American option values from discrete-to continuous-time 
financial models , Mathematical Finance, 4, 289-304. 

[3] L. Andersen (2006): Efficient Simulation of the Heston Stochastic Volatility Model. Preprint available at 
http://www.ressources-actuarielles.net/ 

[4] E. Appolloni, L. Caramellino, A. Zanette (2015): A robust tree method for pricing American options 
with CIR stochastic interest rate. IMA Journal of Management Mathematics, 26, 345-375. 

[5] A. Berman, R. J. Plemmons (1994): Nonnegative matrices in the mathematical sciences, Society for 
Industrial and Applied Mathematics (SIAM), Philadelphia, PA. 

[6] M. Briani, L. Caramellino, A. Zanette (2015): A hybrid approach for the implementation of the 
Heston model. IMA Journal of Management Mathematics, to appear. arXiv: 1307.7178, 


30 



[7] D. Brigo, F. Mercurio (2006): Interest Rate Models-Theory and Practice. Springer, Berlin. 

[8] L. Brugnano, D. Trigiante (1992): Tridiagonal matrices: Invertibility and conditioning, Linear Algebra 
and its Applications, 166, 131-150. 

[9] J.C. Cox, J. Ingersoll, S. Ross (1985): A theory of the term structure of interest rates, Econometrica, 
53, 385-407. 

[10] J. Cox, S.A. Ross, M. Rubinstein (1979): Option pricing: a simplifiled approach. Journal of Financial 
Economics 7, 229-263. 

[11] V. D’Halluin, P.A. Forsyth, G. Labahn (2005): A semi-Lagrangian Approach for American Asian 
options under jump-diffusion, Siam J.Sci.Comp. 27, 315-345. 

[12] S.N. Ethier, T. Kurtz (1986): Markov processes: characterization and convergence. John Wiley & Sons, 
New York. 

[13] T. Haentjens, K.J. in’t Hout (2012): Alternating direction implicit finite difference schemes for the 
Heston-Hull-White partial differential equation. J. Comp. Finan. 16, 83-110. 

[14] J. Hull, A. White A (1994): Numerical procedures for implementing term structure models 1. Journal 
of Derivatives 2(1), 7-16. 

[15] A.L. Grzelak, G.W. Oosterlee (2011): On the Heston model with stochastic interest rates. SIAM J. 
Fin. Math. 2, 255-286. 

[16] A.L. Grzelak, G.W. Oosterlee (2012): On the Cross-currency with stochastic volatility and stochastic 
interest rate. Applied Mathematical Finance 19(1), 1-35 

[17] J.E. Hilliard, A.L. Schwartz, A.L. Tucker (1996): Bivariate binomial pricing with generalized 
interest rate processes. The Journal of Financial Research, XIX-4, 585-602. 

[18] D. Lamberton, G. Pages (1990): Sur I’approximation des reduites. Ann. Inst. H. Poincare, Probab. et 
Statistiques 26, 331-355. 

[19] F.A. Longstaff, E.S. Schwartz (2001): Valuing American options by simulations: a simple least 
squares approach. The Review of Financial Studies, 14, 113-148. 

[20] D.B. Nelson, K. Ramaswamy (1990): Simple binomial processes as diffusion approximations in financial 
models. The Review of Financial Studies, 3, 393-430. 

[21] QuantLib. a free/open-source library for quantitative finance, http://quantlib.org/index.shtml 

[22] M. Vellekoop, H. Nieuwenhuis (2009): A tree-based method to price American Options in the Heston 
Model. The Journal of Computational Finance, 13, 1-21. 

[23] J.Z. Wei (1996): Valuing American equity options with a stochastic interest rate: a note. The Journal of 
Financial Engineering, 2, 195-206. 


31 



